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Abstract. In a recent numerical study [Ng et al., Astrophys. J. 747, 109, 2012], with 
a three-dimensional model of coronal heating using reduced magnetohydrodynamics 
(RMHD), we have obtained scaling results of heating rate versus Lundquist number 
based on a series of runs in which random photospheric motions are imposed for hun- 
dreds to thousands of Alfven time in order to obtain converged statistical values. The 
heating rate found in these simulations saturate to a level that is independent of the 
Lundquist number. This scaling result was also supported by an analysis with the as- 
sumption of the Sweet-Parker scaling of the current sheets, as well as how the width, 
length and number of current sheets scale with Lundquist number In order to test these 
assumptions, we have implemented an automated routine to analyze thousands of cur- 
rent sheets in these simulations and return statistical scalings for these quantities. It is 
found that the Sweet-Parker scaling is justified. However, some discrepancies are also 
found and require further study. 



1. Introduction 

Within the framework of the parker model of coronal heating ('Parker"l972*), a recent 
analysis (Ng & Bhattacharjee 2008) in two dimensions (2D) demonstrated that when 
coherence times (Tc) of the imposed photospheric turbulence are much smaller than 
characteristic resistive time-scales (tr), the Ohmic dissipation scales independently of 
resistivity. While their initial 2D RMHD treatment precluded non-linear effects such 
as instabilities and/or magnetic reconnection, they further invoked a simple analytical 
argument demonstrated that even with these non-linear elfects, which would limit the 
growth of B± (the component of the magnetic field perpendicular to a guide field B^), 
the insensitivity to resistivity is still true for small enough r^.. This latter hypothesis 
was subsequently confirmed (Ng et al.l l2012l) by means of 3D RMHD simulations of 



the Parker model which spanned three orders of magnitude in Lundquist number. 

The scalings derived in these studies depend critically on the assumption that the 
classical 2D steady-state Sweet-Parker scaling for magnetic reconnection holds in 3D 
simulations where extended current sheets form due to random boundary driving. In 
order to empirically substantiate this assumption, and to look into the nature of stochas- 
tically driven magentofluids in more details, we have carried out a systematic analysis 
of current sheets formed in our simulations. 
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Figure 1. 3D iso-surfaces of current sheets. 

We report here a simple algorithm to identify and characterize individual current 
layers. Statistics are accumulated for current sheet parameters. In Section 2, we develop 
and review the motivating analysis. Section 3 describes our simple algorithm. Section 
4 reports our findings and outstanding issues with our analysis. 

2. Coronal Heating Scaling Analysis 

In the parker model, a solar coronal loop is treated as a straight ideal plasma column, 
bounded by two perfectly conducting end-plates representing the photosphere. Initially, 
there is a uniform magnetic field along the z direction. The footpoints of the magnetic 
field on the photosphere are frozen (line-tied), subjected to slow, random motion (p{z = 
0, t) and 0(z = L, t) that deform the magnetic field. 

The footpoint motion is assumed to take place on a time scale much longer than 
the characteristic time for Alfven wave propagation between z - and z = L, so that 
the plasma can be assumed to be in quasi static equilibrium nearly everywhere, if such 
equilibrium exists, during this random evolution. For a given equilibrium, a footpoint 
mapping can be defined by following field lines from one plate to the other. Since the 
plasma is assumed to obey the ideal MHD equations, the magnetic field lines are frozen 
in the plasma and cannot be broken during the twisting process. There f ore, t he foot- 
point mapping must be continuous for smooth footpoint motion. iParkerj (Il972h claims 



that if a sequence of random footpoint motion renders the mapping sufficiently com- 
plicated, there will be no smooth equilibrium for the plasma to relax to, and tangential 
discontinuities (or current sheets) of the magnetic field must develop. Parker treated the 
corona as ideal given that the Lundquist number of the corona is estimated to be of the 
order of 10 Being of ultimately finite resistivity however, it is suggested that ohmic 
dissipation in these current sheets can significantly account for heating coronal plasma. 
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Resolving current sheets at realistic values of Lundquist numbers remains well 
beyond the reach of current computational capabilities with direct numerical MHD 
simulations. Useful scaling studies, however, have been carried out by numerous inves- 
tigators (cf. iNg et al.ll2012l for references). Our current numerical study is motivated 
by scaling analysis for coronal heating that, unlike previous derivations of the coronal 
heating rate, considers the ef fects of random f ootpoint moti on. We quickly summarize 
the main arguments here (cf. iNg & Bhattachar jee 2008 and iNg et al.ll2012 for details). 

In well resolved direct time-dependent, 3D RMHD numerical simulations of the 
Parker model, on average, at any given time, there will be A'^ current sheets with char- 
acteristic width (A) and length (A), which span the length of the simulation box (L). 
The characteristic time over which energy is built-up by random photospheric motions 
and subsequently released is te, and therfore the average heating rate due to ohmic 
dissipation can be written as: 



B 

W ~ rjNAL^ 
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where t] is the plasma resistivity and B\ is the average magnetic field. An expression 
for the average perpendicular magnetic field can be estimated considering that over 
a duration of te, photospheric footpoint motions of average velocity Vp, if assumed 
constant, would deform the guide field B^ and produce B± up to a level of 



B± ~ B, 
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where, for the latter expression, we have solved for te in Equation ([T]) using the Sweet- 



1 /2 

Parker current sheet scaling given hy Aj Is. ~ S I . We define the Lundquist number 
here as S'x = wBj_/t], with w = VpTc being the typical transverse length scale of the 
magnetic field. Together, Equations ([T]) and ^ yield the following expression for the 
average heating rate: 



L'^Bhl\"^ 



(3) 



Evidently, when one extrapolates to the coUisionless coronal limit, the heating rate 
predicted here becomes un-physically large. If we rewrite however the expression for 
the perpendicular magnetic field production considering the turbulent motions in the 
photosphere to have a random walk nature: B^ ~ B^Vp{TcTE)^^^ IL, an average heating 
rate can be estimated as: 



W 



[2 



Equation ^ is manifestly independent of resistivity, a nd holds when 



(4) 



< Te, 



i.e. when the effects of random motion become important. In iNg et all (120121) we have 
studied the transition of the heating rate into this regime with numerical simulations 
spanning three orders of magnitude. In th e regirne whe re Tc > te, we find good agree- 
ment with a previous study (Longcope & Sudanlll994) . with the heating rate scales as 



W oc T]^/^. While in the high Lundquist number regi me wher e Tr < te, we recover 



an T] independent behavior The reader is referred to Ng et all (120121) for an in depth 
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discussion of these results. Here we focus on an ancillary study that addresses two key 
assumptions made in arriving at the scaling relations here: (1) The number of current 
sheets N is essentially independent of Lundquist number; and (2) The Sweet-Parker 
scaling, which derives from a 2D quasi-steady theory, is applicable more generally in 
our 3D line-tied model with driving applied at the boundaries. We assess these two 
assumptions by means of a straightforward algorithm. 



3. Current Sheet Identification and Fitting 

Our MHD simulations employ the reduced M HD equations and are solved using a 
standard algorithm (described in iNg et aL 2012> and ha ve been recently accelerated via 



GPUs with NVidia CUDA as shown in Lin et al.ll2012l) . Our current goal is to identify 



and characterize individual current layers forming in time -dependent 3D simulations of 
a coronal loop driven at the line-tied boundaries. Figure [T]shows iso-surfaces of current 
density for a snapshot of one of our simulation runs (with S j_ ~ 5000). As a starting 
point, we examine only the 2D cross section at z = L/2 of our simulation domain. 

This task is formidable for the following reasons: (1) Given the stochastic nature of 
the imposed photospheric boundary driving, current sheet orientations are random. (2) 
We are dealing with tens of thousands of individual instances of current sheets forming 
during steady state evolution of the Parker model, for which we have data cubes saved 
at a prescribed cadence. (3) We use periodic boundary conditions in which current 
sheets often traverse the edges. (4) In three dimensions, current sheets appear to branch 
out, so a structure appearing as a single current layer in one specific cross-section of 
the loop might appear as several in a different cross-section at a location further along 
the loop, possibly with different properties. Figure [T] attests to each of these issues. 

Our approach consists of two steps. First, an ad-hoc thresholding algorithm iden- 
tifies current sheet candidates by simply taking all pixels in |7| above a pre-defined 
fraction of \J\,„ax and testing for contiguity of the selected regions. This is done in two- 
dimensional cross-sections of the loop simulations, with the algorithm accounting for 
the periodic boundary conditions used by the pseudo-spectral RMHD scheme. By this 
we mean that a current sheet that appears at a border of the simulation box will appear 
at the other border (or at up to 4 edges if it appears at a comer), but will be identified 
only one occurrence. This feature is crucial, considering that we are automating this 
procedure to analyze tens of thousands of simulation cube samples and the likelihood 
of current sheets appearing at domain edges is quite high. Figure |2la) shows a contour 
plot of cunent density for a cross-section of one of our simulations. Current sheet can- 
didates identified by the routine are labeled by green bracketed numbers. The structures 
labeled [1] and [6], for example, appear at edges but are uniquely identified. 

After current sheet candidates are identified they are morphologically examined 
by another automated algorithm, which performs least-square fitting with a bi-variate 
Gaussian. The automated algorithm is implemented using fitting and parameter con- 
straining tool s found in the Package for the Interactive Analysis of Line Emission 
(PINTofALE. lKashvap & D rake 2000). 

Together, these two algorithms yield current sheet orientations with respect to the 
axes (6), number of current sheets present (N), local Jmax, together with CTsmaU and 
ciarge, which scrvc as proxies for current sheet width (A) and current sheet length (A), 
respectively. The three other panels in Figure |2]show one such fit for the current sheet 
labeled [9]. A surface plot shows |/| in the region where current sheet [9] resides in 
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Figure 2. Example of current sheet identification and fitting. 

plot (b). Plot(c) shows the bi-variate Gaussian fit and (d) shows the residual. The pri- 
mary shortfalls of this approach can be summarized as follows: (1) Many current sheets 
are not well approximated by bi-variate Gaussians. Profiles are often asymmetric, and 
the 2-D support of the current sheet structures is often bow- shaped rather than linear. 

(2) Because we are taking only discrete samples in time (full data cubes are saved at 
a pre-determined cadence during simulations runs), the measurements will be biased 
towards current sheet structure that is most long-lived during the lifetime of the sheets. 

(3) In the present analysis, the threshold for structure detection is set at 10% of max- 
imum |7| (as measured in each time step), low enough so that most local maxima are 
included. Unfortunately, this renders the iterative search approach we take, while ro- 
bust, quite computationally expensive. At higher resolutions, this becomes prohibitive, 
and requires re-sampling to lower resolution for reasonable run-times. 

4. Summary of Results & Conclusions 

Our current sheet identification and fitting routines have been applied to high resolution 
data set of Ng et al. (2012), producing a large sample of fits in the range of thousands 
to tens of thousands of current sheets per simulation run. 

In Figure [3l we report weighted average quantities, where we use goodness-of-fit 
from the least square bi-variate Gaussian fitting as the weighting factor. Plots (a) and 
(b) show weighted average values for current sheet widths and lengths as a function 
of Lundquist number S ±. Plot (c) demonstrates good agreement with Sweet-Parker 
scaling. Plot (d) shows the number of current sheets averaged over all post-processed 
time slices. These results provide some empirical support for the assumptions we used 
for our scaling analysis. It is noted however that, even if the Sweet-Parker scaling of 
/}/A~'^^ oc S ^'^^ is recovered, both the current sheet width (Aoc S ^^^) and current sheet 
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Figure 3. Scaling of measured current sheet parameters with Lundquist number. 
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length (A oc ) decrease wi th the inc r ease o f S , faster than the scalings (A x S j_ 
A ~ constant) assumed in both iNg et all ( 2012 ) and lLongcope & Sudani (1994). 

A more detailed analysis of these current sheet fitting results, and an extension of 
this analysis to examine the 3D structure of current layers is now underway. Of partic- 
ular interest here is how the spatial separation of dissipative events in these simulations 
can inform an analysis of flare energy distributions, which typically only consider tem- 
poral variations in event definition (cf. Buchlin et al. 2005; Ng & Lin 20121) . 
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